Load the data and look at data structure.
library(tidyverse)
library(rethinking)
library(tidybayes)
dat = read_csv("clarkia_transplant_data.csv")
dat_nov_germ = list(
nov_germ = dat$nov_germ,
blk = dat$blk,
temperature_diff = dat$temperature_diff_fall,
pop = as.integer(as.factor(dat$pop))
)
table(round(dat_nov_germ$temperature_diff, 2), dat_nov_germ$pop)
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## -1.59 0 0 0 0 236 0 0 0 0 0 0 0 0 0 0
## -1.49 0 0 0 188 0 0 0 0 0 0 0 0 0 0 0
## -1.3 228 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## -1.1 0 0 0 0 0 0 0 226 0 0 0 0 0 0 0
## -0.69 0 0 0 0 0 229 0 0 0 0 0 0 0 0 0
## -0.14 0 0 0 0 0 0 0 0 0 0 0 0 0 230 0
## -0.04 0 0 234 0 0 0 0 0 0 0 0 0 0 0 0
## 0.34 0 0 0 0 0 0 0 0 0 0 176 0 0 0 0
## 0.68 0 0 0 0 0 0 0 0 0 0 0 0 0 0 235
## 0.7 0 236 0 0 0 0 0 0 0 0 0 0 0 0 0
## 1.19 0 0 0 0 0 0 200 0 0 0 0 0 0 0 0
## 1.54 0 0 0 0 0 0 0 0 0 0 0 240 0 0 0
## 1.6 0 0 0 0 0 0 0 0 0 0 0 0 128 0 0
## 1.67 0 0 0 0 0 0 0 0 232 0 0 0 0 0 0
## 1.94 0 0 0 0 0 0 0 0 0 232 0 0 0 0 0
table(dat_nov_germ$blk, dat_nov_germ$pop)
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 1 56 59 59 47 59 58 50 55 58 57 44 60 32 58 58
## 2 56 59 57 47 59 56 50 57 58 57 44 60 32 56 59
## 3 58 59 59 47 59 57 50 57 58 59 44 60 32 58 59
## 4 58 59 59 47 59 58 50 57 58 59 44 60 32 58 59
Fit model:
# Model with temperature only
m0 = ulam(
alist(
nov_germ ~ dbinom(1, p),
logit(p) <- a + bT*temperature_diff,
a ~ dnorm(0, 1.5),
bT ~ dnorm(0, 0.5)
), data = dat_nov_germ, cores = 2, chains = 2, iter = 2000, log_lik = TRUE)
precis(m0, depth = 2)
## mean sd 5.5% 94.5% n_eff Rhat
## a 0.2047734 0.03642652 0.1456934 0.2596579 1216.346 1.001950
## bT -0.2127429 0.02973011 -0.2602620 -0.1661672 1057.771 1.002984
Check on chains:
traceplot(m0)
prior = extract.prior(m0)
##
## SAMPLING FOR MODEL 'c0628fc6583a8e27a955a78f4cdc6094' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 0.000493 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.93 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.23158 seconds (Warm-up)
## Chain 1: 0.957521 seconds (Sampling)
## Chain 1: 2.1891 seconds (Total)
## Chain 1:
mu = link(m0, post = prior, data = list(temperature_diff = seq(-1.5, 2, by = 0.5)))
mu_plot = data.frame(mu) %>%
mutate(row = row_number()) %>%
pivot_longer(cols = X1:X8, names_to = "key", values_to = "value") %>%
mutate(temp_diff = (as.numeric(str_remove(key, "X")) - 4)/2)
ggplot(mu_plot, aes(x = temp_diff, y = value, group = row)) +
geom_line(alpha = 0.2)
Seems ok, allows for many relationships between temperature difference and germination.